options(scipen = 999, digits = 4)
knitr::opts_chunk$set(comment = "#")
r = getOption("repos")
r["CRAN"] = "http://cran.us.r-project.org"
options(repos = r)
## load required packages
ipak <- function (pkg) {
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if (length(new.pkg))
install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c("tidyverse", "Hmisc", "lattice", "multcomp", "lsmeans", "schoRsch", "influence.ME", "lme4", "effects", "lmerTest", "cowplot", "irr", "simr", "plyr", "dplyr","patchwork", "wesanderson", "MuMIn", "devtools", "dplyr", "ggResidpanel", "HLMdiag", "mixed", "sjPlot")
ipak(packages)
## Warning: package 'mixed' is not available (for R version 4.0.2)
## Warning: replacing previous import 'vctrs::data_frame' by 'tibble::data_frame'
## when loading 'dplyr'
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'mixed'
## tidyverse Hmisc lattice multcomp lsmeans schoRsch
## TRUE TRUE TRUE TRUE TRUE TRUE
## influence.ME lme4 effects lmerTest cowplot irr
## TRUE TRUE TRUE TRUE TRUE TRUE
## simr plyr dplyr patchwork wesanderson MuMIn
## TRUE TRUE TRUE TRUE TRUE TRUE
## devtools dplyr ggResidpanel HLMdiag mixed sjPlot
## TRUE TRUE TRUE TRUE FALSE TRUE
# set summed contrasts
options(contrasts = c("contr.sum", "contr.poly"))
# fix bs with dplyr
detach("package:dplyr", unload = TRUE)
## Warning: 'dplyr' namespace cannot be unloaded:
## namespace 'dplyr' is imported by 'tidyr', 'sjPlot', 'plotly', 'broom', 'janitor', 'dbplyr', 'sjmisc', 'sjstats', 'qqplotr', 'HLMdiag' so cannot be unloaded
library("dplyr")
ind.data.pre <- read.csv("./individual_data.csv") %>%
mutate(experiment=str_sub(condition,1,1)) %>%
mutate_if(is.character,as.factor) %>%
mutate(subj = as.factor(subj))%>%
mutate_at(.vars="condition", .funs = tolower)
ma.data <- read.csv("./ma_data.csv") %>%
rename(paper = study_ID,
condition = expt_condition,
experiment = expt_num,
task = looking_task) %>%
mutate_if(is.character,as.factor) %>%
mutate_at(.vars="condition", .funs = tolower)
study.design <- ma.data %>%
select(paper, experiment, condition, task, training_yesno, action_consequence, actor_hand, agent_efficient_fam, object_diff_size_huge, action_causal, action_consequence, location_object_goal_ambiguous, bothobjects_present_visible_fam, agent, PI_group) %>%
mutate_if(is.character,as.factor) %>%
mutate(experiment = as.factor(experiment)) %>%
filter(paper != "sommerville2005",
task != "causes")
# merge study design with individual looks from babies
ind.data <- left_join(ind.data.pre, study.design,
by=c("paper", "experiment", "condition")) %>%
mutate(sex = tolower(str_sub(sex,1,1)),
look_pref = unexp_look - exp_look) %>%
mutate(paper = str_replace_all(paper, "unpublished", "unpub")) %>%
mutate(task = str_replace_all(task, "efficiency", "constraints")) %>%
mutate_if(is.character,as.factor)
str(ind.data)
# 'data.frame': 626 obs. of 20 variables:
# $ subj : Factor w/ 626 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
# $ condition : Factor w/ 25 levels "1_active","1_control",..: 12 12 12 12 12 12 12 12 12 12 ...
# $ ageday : num 103.3 86.9 105.3 129.8 107.3 ...
# $ exp_look : num 13.1 18.4 15.1 7.1 5.8 9 18.9 6 31 42.9 ...
# $ unexp_look : num 60 60 46.3 13.6 12 60 8 21 60 42.7 ...
# $ sex : Factor w/ 2 levels "f","m": 2 2 2 1 1 1 2 2 2 1 ...
# $ paper : Factor w/ 8 levels "choi_unpub","choi2018",..: 6 6 6 6 6 6 6 6 6 6 ...
# $ experiment : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
# $ task : Factor w/ 2 levels "constraints",..: 2 2 2 2 2 2 2 2 2 2 ...
# $ training_yesno : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
# $ action_consequence : Factor w/ 3 levels "location_change",..: 2 2 2 2 2 2 2 2 2 2 ...
# $ actor_hand : Factor w/ 3 levels "bare","gloved",..: NA NA NA NA NA NA NA NA NA NA ...
# $ agent_efficient_fam : Factor w/ 2 levels "no","yes": NA NA NA NA NA NA NA NA NA NA ...
# $ object_diff_size_huge : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
# $ action_causal : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
# $ location_object_goal_ambiguous : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
# $ bothobjects_present_visible_fam: Factor w/ 3 levels "","no","yes": 3 3 3 3 3 3 3 3 3 3 ...
# $ agent : Factor w/ 3 levels "hand","object",..: 2 2 2 2 2 2 2 2 2 2 ...
# $ PI_group : Factor w/ 5 levels "desrochers","luo",..: 2 2 2 2 2 2 2 2 2 2 ...
# $ look_pref : num 46.9 41.6 31.2 6.5 6.2 ...
# View(ind.data)
# function that returns column of standardized betas from lmer model
gen.beta <- function(model) {
df <- data.frame(fixef(model))
names(df) <- c("beta")
return(df)
}
# function that computes CIs and returns them in df
gen.ci <- function(model) {
df <- data.frame(confint(model))
names(df) <- c("lower", "upper")
return(df)
}
# function that converts model summary (lmer) to df
gen.m <- function(model) {
df <- data.frame(coef(summary(model)))
names(df) <- c("b", "se", "df", "t", "p")
return(df)
}
# function that converts model summary (lm) to df
gen.lm <- function(model) {
df <- data.frame(coef(summary(model)))
names(df) <- c("b", "se", "t", "p")
return(df)
}
# function that returns age info and number of females in a dataset
ages <- function(longdata) {
longdata %>% summarize(mean = mean(ageday), min=range(ageday)[1], max=range(ageday)[2], f=sum(sex=="f")/2)
}
# function that returns formatted result from lme4/lmerTest table
report <- function(table, index, places, tails, flip) {
if (tails == "1") {
p <- round(table$p[index], places)/2
howmanytails <- "one-tailed"
} else {
p <- round(table$p[index], places)
howmanytails <- "two-tailed"
}
if (p < .001) {
p <- "<.001"
} else {
p <- paste("=", round(p, places), sep = "")
}
if (missing(flip)) {
result <- paste("[", round(table$lower[index], places), ",", round(table$upper[index], places), "], ß=", round(table$beta[index], places), ", B=", round(table$b[index],places), ", SE=", round(table$se[index],places), ", p", p, ", ", howmanytails, sep = "")
} else {
result <- paste("[", -round(table$upper[index], places), ",", -round(table$lower[index], places), "], ß=", -round(table$beta[index], places), ", B=", -round(table$b[index],places), ", SE=", round(table$se[index],places), ", p", p, ", ", howmanytails, sep = "")
}
return(result)
}
describe <- function(dataset){
summary <- dataset %>%
summarise(unexp_avg = mean(unexp_look),
exp_avg = mean(exp_look),
lookdiff_avg = mean(look_pref),
n = n())
paste(
#"M_unexp = ", summary$unexp_avg, "s ",
#"M_exp = ", summary$exp_avg, "s ",
"Mean_diff = ", round(summary$lookdiff_avg,3), "seconds"
)
}
## Retrieved from : http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/#error-bars-for-within-subjects-variables
## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
## data: a data frame.
## measurevar: the name of a column that contains the variable to be summariezed
## groupvars: a vector containing names of columns that contain grouping variables
## na.rm: a boolean that indicates whether to ignore NA's
## conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
conf.interval=.95, .drop=TRUE) {
library(plyr)
# New version of length which can handle NA's: if na.rm==T, don't count them
length2 <- function (x, na.rm=FALSE) {
if (na.rm) sum(!is.na(x))
else length(x)
}
# This does the summary. For each group's data frame, return a vector with
# N, mean, and sd
datac <- ddply(data, groupvars, .drop=.drop,
.fun = function(xx, col) {
c(N = length2(xx[[col]], na.rm=na.rm),
mean = mean (xx[[col]], na.rm=na.rm),
sd = sd (xx[[col]], na.rm=na.rm)
)
},
measurevar
)
# Rename the "mean" column
datac <- plyr::rename(datac, c("mean" = measurevar))
datac$se <- datac$sd / sqrt(datac$N) # Calculate standard error of the mean
# Confidence interval multiplier for standard error
# Calculate t-statistic for confidence interval:
# e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
ciMult <- qt(conf.interval/2 + .5, datac$N-1)
datac$ci <- datac$se * ciMult
return(datac)
}
## Norms the data within specified groups in a data frame; it normalizes each
## subject (identified by idvar) so that they have the same mean, within each group
## specified by betweenvars.
## data: a data frame.
## idvar: the name of a column that identifies each subject (or matched subjects)
## measurevar: the name of a column that contains the variable to be summariezed
## betweenvars: a vector containing names of columns that are between-subjects variables
## na.rm: a boolean that indicates whether to ignore NA's
normDataWithin <- function(data=NULL, idvar, measurevar, betweenvars=NULL,
na.rm=FALSE, .drop=TRUE) {
library(plyr)
# Measure var on left, idvar + between vars on right of formula.
data.subjMean <- ddply(data, c(idvar, betweenvars), .drop=.drop,
.fun = function(xx, col, na.rm) {
c(subjMean = mean(xx[,col], na.rm=na.rm))
},
measurevar,
na.rm
)
# Put the subject means with original data
data <- merge(data, data.subjMean)
# Get the normalized data in a new column
measureNormedVar <- paste(measurevar, "_norm", sep="")
data[,measureNormedVar] <- data[,measurevar] - data[,"subjMean"] +
mean(data[,measurevar], na.rm=na.rm)
# Remove this subject mean column
data$subjMean <- NULL
return(data)
}
## Summarizes data, handling within-subjects variables by removing inter-subject variability.
## It will still work if there are no within-S variables.
## Gives count, un-normed mean, normed mean (with same between-group mean),
## standard deviation, standard error of the mean, and confidence interval.
## If there are within-subject variables, calculate adjusted values using method from Morey (2008).
## data: a data frame.
## measurevar: the name of a column that contains the variable to be summariezed
## betweenvars: a vector containing names of columns that are between-subjects variables
## withinvars: a vector containing names of columns that are within-subjects variables
## idvar: the name of a column that identifies each subject (or matched subjects)
## na.rm: a boolean that indicates whether to ignore NA's
## conf.interval: the percent range of the confidence interval (default is 95%)
summarySEwithin <- function(data=NULL, measurevar, betweenvars=NULL, withinvars=NULL,
idvar=NULL, na.rm=FALSE, conf.interval=.95, .drop=TRUE) {
# Ensure that the betweenvars and withinvars are factors
factorvars <- vapply(data[, c(betweenvars, withinvars), drop=FALSE],
FUN=is.factor, FUN.VALUE=logical(1))
if (!all(factorvars)) {
nonfactorvars <- names(factorvars)[!factorvars]
message("Automatically converting the following non-factors to factors: ",
paste(nonfactorvars, collapse = ", "))
data[nonfactorvars] <- lapply(data[nonfactorvars], factor)
}
# Get the means from the un-normed data
datac <- summarySE(data, measurevar, groupvars=c(betweenvars, withinvars),
na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)
# Drop all the unused columns (these will be calculated with normed data)
datac$sd <- NULL
datac$se <- NULL
datac$ci <- NULL
# Norm each subject's data
ndata <- normDataWithin(data, idvar, measurevar, betweenvars, na.rm, .drop=.drop)
# This is the name of the new column
measurevar_n <- paste(measurevar, "_norm", sep="")
# Collapse the normed data - now we can treat between and within vars the same
ndatac <- summarySE(ndata, measurevar_n, groupvars=c(betweenvars, withinvars),
na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)
# Apply correction from Morey (2008) to the standard error and confidence interval
# Get the product of the number of conditions of within-S variables
nWithinGroups <- prod(vapply(ndatac[,withinvars, drop=FALSE], FUN=nlevels,
FUN.VALUE=numeric(1)))
correctionFactor <- sqrt( nWithinGroups / (nWithinGroups-1) )
# Apply the correction factor
ndatac$sd <- ndatac$sd * correctionFactor
ndatac$se <- ndatac$se * correctionFactor
ndatac$ci <- ndatac$ci * correctionFactor
# Combine the un-normed means with the normed results
merge(datac, ndatac)
}
ind.data.long <- ind.data %>%
rename(expected = exp_look,
unexpected = unexp_look) %>%
pivot_longer(cols = c(expected, unexpected),
names_to = "trial_type",
values_to = "looking_time")
ind.data.summary <- summarySEwithin(data = ind.data.long, measurevar = "looking_time", betweenvars = c("task", "paper", "experiment", "condition"), withinvars = "trial_type")
# Automatically converting the following non-factors to factors: trial_type
n_infants <- ind.data.summary %>%
filter(trial_type == "expected") %>% # just take one row per study
summarise(sum(N))
n_conditions <- ind.data.summary %>%
filter(trial_type == "expected") %>% # just take one row per study
nrow()
n_papers <- ind.data.summary %>%
filter(trial_type == "expected") %>% # just take one row per study
group_by(paper) %>%
count(paper) %>%
nrow()
We searched for all journal papers, theses, and conference proceedings that reported looking time data from typically developing infants between 2 and 4 months of age in a task structured like the goals or constraints task. We also emailed two listservs (Cognitive Development Society, and Infant Studies) to request more datasets. This resulted in a final list of 11 papers and 35 conditions. We contacted the authors and asked them to send us the original datasets from this past published and unpublished work. We were able to gather original datasets from 8 papers, 29 conditions, and 626 infants (age range 72-137. A team of researchers were then randomly assigned to look through relevant papers including supplemental materials (with SL double-checking every entry and resolving disagreements). When exact values were not provided, or when we came across other ambiguities (e.g. numbers reported in the paper differ from the numbers calculated using the raw data), the team contacted authors to try and address, and also used the tool WebplotDigitizer (https://automeris.io/WebPlotDigitizer/) to extract estimated values from figures. If there were discrepancies between the paper, figures and/or raw data, we prioritized the values from the raw data if available, then author correspondence, then paper, then estimates from figures, in that order. For individual datasets, authors were asked to provide their data and a codebook, and were asked for permission to share their stimuli and data publicly on OSF. Of 8 papers (29, conditions), data from 25 conditions and stimuli from 13 conditions (either actual study videos, or example stimuli) are publicly available at https://osf.io/zwncg/.
table.goals <- ind.data %>%
filter(task=="goals") %>%
group_by(paper, condition, training_yesno, action_causal, action_consequence, location_object_goal_ambiguous, bothobjects_present_visible_fam, agent) %>%
summarise(n = n(), agemin = floor(range(ageday)[1]), agemax = floor(range(ageday)[2])) %>%
select(paper, condition, n, agemin, agemax, training_yesno, action_causal, action_consequence, location_object_goal_ambiguous, bothobjects_present_visible_fam, agent)
# `summarise()` regrouping output by 'paper', 'condition', 'training_yesno', 'action_causal', 'action_consequence', 'location_object_goal_ambiguous', 'bothobjects_present_visible_fam' (override with `.groups` argument)
table.constraints <- ind.data %>%
filter(task=="constraints") %>%
group_by(paper, condition, training_yesno, action_causal, action_consequence, location_object_goal_ambiguous, agent_efficient_fam, actor_hand) %>%
summarise(n = n(), agemin = floor(range(ageday)[1]), agemax = floor(range(ageday)[2])) %>%
select(paper, condition, n, agemin, agemax, training_yesno, action_causal, action_consequence, agent_efficient_fam, actor_hand)
# `summarise()` regrouping output by 'paper', 'condition', 'training_yesno', 'action_causal', 'action_consequence', 'location_object_goal_ambiguous', 'agent_efficient_fam' (override with `.groups` argument)
# Adding missing grouping variables: `location_object_goal_ambiguous`
knitr::kable(table.goals)
| paper | condition | n | agemin | agemax | training_yesno | action_causal | action_consequence | location_object_goal_ambiguous | bothobjects_present_visible_fam | agent |
|---|---|---|---|---|---|---|---|---|---|---|
| choi_unpub | 1_equidistant | 24 | 76 | 137 | no | no | none | yes | yes | person |
| choi_unpub | 1_fartarget | 24 | 72 | 130 | no | no | none | yes | no | person |
| choi_unpub | 1_neartarget | 24 | 74 | 127 | no | no | none | yes | yes | person |
| choi2018 | 1_hidden | 16 | 75 | 121 | no | no | none | yes | yes | person |
| choi2018 | 1_oneobject | 16 | 75 | 127 | no | no | none | yes | yes | person |
| choi2018 | 1_twoobject | 16 | 77 | 130 | no | no | none | yes | yes | person |
| gerson2014a | 1_active | 24 | 91 | 121 | yes | no | none | yes | yes | hand |
| gerson2014a | 1_control | 24 | 91 | 120 | no | no | none | yes | yes | hand |
| gerson2014a | 1_observation | 24 | 91 | 121 | no | no | none | yes | yes | hand |
| gerson2014b | 1_active | 30 | 91 | 125 | yes | no | none | yes | yes | hand |
| gerson2014b | 1_observation | 30 | 92 | 125 | no | no | none | yes | yes | hand |
| gerson2014b | 2_generalization | 30 | 94 | 125 | yes | no | none | yes | yes | hand |
| luo2011 | 1_oneobject | 12 | 76 | 124 | no | no | none | yes | no | object |
| luo2011 | 1_twoobject | 12 | 79 | 129 | no | no | none | yes | yes | object |
| luo2011 | 2_differentpositions | 12 | 81 | 118 | no | no | none | no | no | object |
| woo_unpub | 1_statechange_objectswap | 20 | 93 | 120 | no | yes | state_change | yes | yes | person |
| woo_unpub | 2_disambiguatingobjectgoal | 24 | 91 | 121 | no | yes | state_change | no | yes | person |
knitr::kable(table.constraints)
| location_object_goal_ambiguous | paper | condition | n | agemin | agemax | training_yesno | action_causal | action_consequence | agent_efficient_fam | actor_hand |
|---|---|---|---|---|---|---|---|---|---|---|
| NA | liu2019 | 1_pickupglove | 20 | 92 | 122 | no | yes | location_change | yes | gloved |
| NA | liu2019 | 2_pickupbarehand | 20 | 93 | 120 | no | yes | location_change | yes | bare |
| NA | liu2019 | 3_statechangebarrier | 20 | 91 | 122 | no | yes | state_change | yes | gloved |
| NA | liu2019 | 3_statechangenobarrier | 20 | 91 | 122 | no | yes | state_change | no | gloved |
| NA | liu2019 | 4_statechangenotcausal | 20 | 93 | 121 | no | no | state_change | yes | gloved |
| NA | liu2019 | 5_statechangecausal | 26 | 92 | 121 | no | yes | state_change | yes | gloved |
| NA | liu2019 | 5_statechangenotcausal | 26 | 93 | 120 | no | no | state_change | yes | gloved |
| NA | skerry2013 | 1_effectiveactiontraining | 20 | 93 | 121 | yes | yes | location_change | yes | mittened |
| NA | skerry2013 | 2_ineffectiveactiontraining | 20 | 93 | 122 | no | yes | location_change | yes | mittened |
| NA | skerry2013 | 3_notraining | 20 | 91 | 121 | no | yes | location_change | yes | mittened |
| NA | skerry2013 | 4_constrainedactionhabituation | 26 | 93 | 120 | yes | yes | location_change | yes | mittened |
| NA | skerry2013 | 5_unconstrainedactionhabituation | 26 | 93 | 122 | yes | yes | location_change | no | mittened |
str(ind.data)
# 'data.frame': 626 obs. of 20 variables:
# $ subj : Factor w/ 626 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
# $ condition : Factor w/ 25 levels "1_active","1_control",..: 12 12 12 12 12 12 12 12 12 12 ...
# $ ageday : num 103.3 86.9 105.3 129.8 107.3 ...
# $ exp_look : num 13.1 18.4 15.1 7.1 5.8 9 18.9 6 31 42.9 ...
# $ unexp_look : num 60 60 46.3 13.6 12 60 8 21 60 42.7 ...
# $ sex : Factor w/ 2 levels "f","m": 2 2 2 1 1 1 2 2 2 1 ...
# $ paper : Factor w/ 8 levels "choi_unpub","choi2018",..: 6 6 6 6 6 6 6 6 6 6 ...
# $ experiment : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
# $ task : Factor w/ 2 levels "constraints",..: 2 2 2 2 2 2 2 2 2 2 ...
# $ training_yesno : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
# $ action_consequence : Factor w/ 3 levels "location_change",..: 2 2 2 2 2 2 2 2 2 2 ...
# $ actor_hand : Factor w/ 3 levels "bare","gloved",..: NA NA NA NA NA NA NA NA NA NA ...
# $ agent_efficient_fam : Factor w/ 2 levels "no","yes": NA NA NA NA NA NA NA NA NA NA ...
# $ object_diff_size_huge : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
# $ action_causal : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
# $ location_object_goal_ambiguous : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
# $ bothobjects_present_visible_fam: Factor w/ 3 levels "","no","yes": 3 3 3 3 3 3 3 3 3 3 ...
# $ agent : Factor w/ 3 levels "hand","object",..: 2 2 2 2 2 2 2 2 2 2 ...
# $ PI_group : Factor w/ 5 levels "desrochers","luo",..: 2 2 2 2 2 2 2 2 2 2 ...
# $ look_pref : num 46.9 41.6 31.2 6.5 6.2 ...
ggplot(ind.data, aes(x=look_pref)) +
geom_histogram()+
facet_grid(~paper)+
theme_cowplot(12)
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
plot1 <- ggplot(ind.data.long,
aes(trial_type, looking_time, colour=paper))+
theme_cowplot(12)+
theme(legend.title = element_blank())+
geom_boxplot(outlier.alpha = 0.2)+
geom_point(alpha=0.2)+
geom_line(aes(group=subj), alpha=0.2)+
stat_summary(geom="point", fun="mean", colour="black")+
# geom_errorbar(data = ind.data.summary, colour="red", position = position_dodge(width = 5), width = 0, aes(ymin=looking_time-se, ymax=looking_time+se)) +
stat_summary(geom="errorbar", fun.data="mean_se", width=0.2, colour="black")+
ylab("Looking Time (s)")+
xlab("Condition")+
facet_wrap(~task+paper+condition, nrow=4, labeller=label_wrap_gen())+
theme(axis.text.x = element_text(angle = 90, hjust=0.95, vjust =0.2),
legend.position="bottom")
plot1
plot2.goals <- ggplot(ind.data %>% filter(task == "goals"),
aes(condition, look_pref, fill=paper)) +
geom_boxplot(outlier.alpha = 0.2, position = position_dodge2(preserve = "single"))+
geom_hline(yintercept = 0)+
geom_point(alpha=0.2)+
stat_summary(geom="point", fun="mean", colour="black")+
stat_summary(geom="errorbar", fun.data="mean_se", width=0.2, colour="black")+
theme_cowplot(18)+
facet_grid(~paper, scales = "free_x",space = "free", labeller = label_wrap_gen())+
ylab("Looking preference (s) \n <- Expected ------- Unexpected ->")+
theme(axis.text.x = element_text(angle = 90, hjust=0.95, vjust =0.2), legend.position = "none")+
scale_fill_brewer(palette = "Dark2")
plot2.constraints <- ggplot(ind.data %>% filter(task == "constraints"),
aes(condition, look_pref, fill=paper)) +
geom_boxplot(outlier.alpha = 0.2, position = position_dodge2(preserve = "single"))+
geom_hline(yintercept = 0)+
geom_point(alpha=0.2)+
stat_summary(geom="point", fun="mean", colour="black")+
stat_summary(geom="errorbar", fun.data="mean_se", width=0.2, colour="black")+
theme_cowplot(18)+
facet_grid(~paper,scales = "free_x",space = "free", labeller = label_wrap_gen())+
ylab("")+
xlab("")+
theme(axis.text.x = element_text(angle = 90, hjust=0.95, vjust =0.2), legend.position = "none")+
scale_fill_manual(values = wes_palette("Royal2"))
plot2.goals
plot2.constraints
(plot2.goals | plot2.constraints) + plot_layout(widths = c(2, 1)) + plot_annotation(tag_levels = 'A')
(plot3 <- ggplot(ind.data,
aes(ageday, look_pref, colour=paper)) +
# geom_boxplot(outlier.alpha = 0.2, position = position_dodge2(preserve = "single"))+
geom_hline(yintercept = 0)+
geom_point()+
geom_smooth(method="lm")+
# stat_summary(geom="point", fun="mean", colour="black")+
# stat_summary(geom="errorbar", fun.data="mean_se", width=0.2, colour="black")+
theme_cowplot(12)+
facet_grid(~task+paper)+
ylab("Looking preference (s) \n <- Expected ------- Unexpected ->")+
theme(axis.text.x = element_text(angle = 90, hjust=0.95, vjust =0.2)))
# `geom_smooth()` using formula 'y ~ x'
constraints <- ind.data %>% filter(task =="constraints")
constraints$actor_hand
# [1] mittened mittened mittened mittened mittened mittened mittened mittened
# [9] mittened mittened mittened mittened mittened mittened mittened mittened
# [17] mittened mittened mittened mittened mittened mittened mittened mittened
# [25] mittened mittened mittened mittened mittened mittened mittened mittened
# [33] mittened mittened mittened mittened mittened mittened mittened mittened
# [41] mittened mittened mittened mittened mittened mittened mittened mittened
# [49] mittened mittened mittened mittened mittened mittened mittened mittened
# [57] mittened mittened mittened mittened mittened mittened mittened mittened
# [65] mittened mittened mittened mittened mittened mittened mittened mittened
# [73] mittened mittened mittened mittened mittened mittened mittened mittened
# [81] mittened mittened mittened mittened mittened mittened mittened mittened
# [89] mittened mittened mittened mittened mittened mittened mittened mittened
# [97] mittened mittened mittened mittened mittened mittened mittened mittened
# [105] mittened mittened mittened mittened mittened mittened mittened mittened
# [113] gloved gloved gloved gloved gloved gloved gloved gloved
# [121] gloved gloved gloved gloved gloved gloved gloved gloved
# [129] gloved gloved gloved gloved gloved gloved gloved gloved
# [137] gloved gloved gloved gloved gloved gloved gloved gloved
# [145] gloved gloved gloved gloved gloved gloved gloved gloved
# [153] gloved gloved gloved gloved gloved gloved gloved gloved
# [161] gloved gloved gloved gloved gloved gloved gloved gloved
# [169] gloved gloved gloved gloved gloved gloved gloved gloved
# [177] gloved gloved gloved gloved gloved gloved gloved gloved
# [185] gloved gloved gloved gloved gloved gloved gloved gloved
# [193] gloved gloved gloved gloved gloved gloved gloved gloved
# [201] gloved gloved gloved gloved gloved gloved gloved gloved
# [209] gloved gloved gloved gloved gloved gloved gloved gloved
# [217] gloved gloved gloved gloved gloved gloved gloved gloved
# [225] gloved gloved gloved gloved gloved gloved gloved gloved
# [233] gloved gloved gloved gloved gloved gloved gloved gloved
# [241] gloved gloved gloved gloved bare bare bare bare
# [249] bare bare bare bare bare bare bare bare
# [257] bare bare bare bare bare bare bare bare
# Levels: bare gloved mittened
constraints$action_consequence <- relevel(constraints$action_consequence, ref = "none")
constraints$agent_efficient_fam <- relevel(constraints$agent_efficient_fam, ref = "no")
constraints$training_yesno <- relevel(constraints$training_yesno, ref = "no")
constraints$action_causal <- relevel(constraints$action_causal, ref = "no")
constraints.baseline.data <- constraints %>%
filter(condition %in% c("3_notraining",
"1_pickupglove",
"2_pickupbarehand"))
constraints.b1 <- lmer(data = constraints.baseline.data,
formula = look_pref ~ 1 + ageday + (1|condition))
summary(constraints.b1)
# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula: look_pref ~ 1 + ageday + (1 | condition)
# Data: constraints.baseline.data
#
# REML criterion at convergence: 366.4
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -2.7189 -0.5297 -0.0435 0.7252 2.2686
#
# Random effects:
# Groups Name Variance Std.Dev.
# condition (Intercept) 2.03 1.42
# Residual 25.37 5.04
# Number of obs: 60, groups: condition, 3
#
# Fixed effects:
# Estimate Std. Error df t value Pr(>|t|)
# (Intercept) -2.3930 9.0407 57.2416 -0.26 0.79
# ageday 0.0258 0.0832 56.0050 0.31 0.76
#
# Correlation of Fixed Effects:
# (Intr)
# ageday -0.993
constraints.b1.table <- sjPlot:: tab_model(constraints.b1,
show.std =TRUE,
show.stat=TRUE,
show.df=TRUE)
constraints.b1.beta <- summary(effectsize::standardize(constraints.b1))
constraints.baseline <- cbind(
gen.beta(effectsize::standardize(constraints.b1)),
gen.m(constraints.b1),
gen.ci(constraints.b1)[3:4,]
)
# Computing profile confidence intervals ...
constraints.baseline
# beta b se df t p
# (Intercept) -0.00000000000000002645 -2.39296 9.04071 57.24 -0.2647 0.7922
# ageday 0.03966247939082377660 0.02582 0.08316 56.01 0.3105 0.7574
# lower upper
# (Intercept) -20.1694 15.4724
# ageday -0.1389 0.1897
cooks1 <- cooks.distance(constraints.b1, group = "subj")
# Warning in cooks.distance.lmerMod(constraints.b1, group = "subj"): group is not
# a valid argument for this function. As of version 0.4.0, group has been replaced
# by level.
dotplot_diag(x = cooks1, cutoff = "internal", name = "cooks.distance", index=constraints.baseline.data$subj) + ylab("Cook's distance") + xlab("subjID")
excluded.subs <- c("90", "553", "86")
constraints.b1.cooks <- lmer(data = constraints.baseline.data %>%
filter(!subj %in% excluded.subs),
formula = look_pref ~ 1 + ageday + (1|condition))
summary(constraints.b1.cooks)
# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula: look_pref ~ 1 + ageday + (1 | condition)
# Data: constraints.baseline.data %>% filter(!subj %in% excluded.subs)
#
# REML criterion at convergence: 338.6
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -2.5544 -0.5380 -0.0717 0.7259 2.4971
#
# Random effects:
# Groups Name Variance Std.Dev.
# condition (Intercept) 1.14 1.07
# Residual 21.61 4.65
# Number of obs: 57, groups: condition, 3
#
# Fixed effects:
# Estimate Std. Error df t value Pr(>|t|)
# (Intercept) -9.5734 9.0731 54.2948 -1.06 0.30
# ageday 0.0926 0.0838 53.5855 1.11 0.27
#
# Correlation of Fixed Effects:
# (Intr)
# ageday -0.995
plot(allEffects(constraints.b1.cooks))
constraints.baseline.cooks <- cbind(
gen.beta(effectsize::standardize(constraints.b1.cooks)),
gen.m(constraints.b1.cooks),
gen.ci(constraints.b1.cooks)[3:4,]
)
# Computing profile confidence intervals ...
For standard versions of the constraints task, we found no differential looking between the expected and unexpected events, Mean_diff = 0.395 seconds, [-0.139,0.19], ß=0.04, B=0.026, SE=0.083, p=0.757, two-tailed. This result held when we removed 3 influential observations, identified using Cook’s distance, [-0.07,0.263], ß=0.145, B=0.093, SE=0.084, p=0.274, two-tailed.
constraints.m1 <- lmer(data = constraints,
formula = look_pref ~ training_yesno + action_causal + action_consequence + actor_hand + agent_efficient_fam + ageday + (1|condition) + (1|experiment) + (1|paper))
# boundary (singular) fit: see ?isSingular
# Warning in as_lmerModLT(model, devfun): Model may not have converged with 1
# eigenvalue close to zero: -3.1e-11
summary(constraints.m1)
# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula: look_pref ~ training_yesno + action_causal + action_consequence +
# actor_hand + agent_efficient_fam + ageday + (1 | condition) +
# (1 | experiment) + (1 | paper)
# Data: constraints
#
# REML criterion at convergence: 1682
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -3.729 -0.510 -0.059 0.561 3.866
#
# Random effects:
# Groups Name Variance Std.Dev.
# condition (Intercept) 0.000 0.000
# experiment (Intercept) 0.000 0.000
# paper (Intercept) 0.121 0.347
# Residual 35.324 5.943
# Number of obs: 264, groups: condition, 12; experiment, 5; paper, 2
#
# Fixed effects:
# Estimate Std. Error df t value Pr(>|t|)
# (Intercept) -2.8929 4.6171 256.0000 -0.63 0.53151
# training_yesno1 -2.1859 0.6174 256.0000 -3.54 0.00047 ***
# action_causal1 -2.3185 0.5944 256.0000 -3.90 0.00012 ***
# action_consequence1 -1.0693 0.7763 256.0000 -1.38 0.16959
# actor_hand1 1.3018 1.0518 256.0000 1.24 0.21698
# actor_hand2 0.8084 1.0518 256.0000 0.77 0.44282
# agent_efficient_fam1 -1.7116 0.5377 256.0000 -3.18 0.00164 **
# ageday 0.0231 0.0419 256.0000 0.55 0.58139
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Correlation of Fixed Effects:
# (Intr) trnn_1 actn_cs1 actn_cn1 actr_1 actr_2 agn__1
# tranng_ysn1 -0.062
# action_csl1 0.143 0.085
# actn_cnsqn1 0.010 0.065 0.349
# actor_hand1 0.124 -0.227 0.001 -0.360
# actor_hand2 -0.060 -0.227 0.000 0.721 -0.597
# agnt_ffcn_1 0.092 0.314 0.275 0.210 0.000 0.000
# ageday -0.982 0.017 -0.053 -0.036 -0.010 -0.008 -0.018
# convergence code: 0
# boundary (singular) fit: see ?isSingular
resid_panel(constraints.m1, plots = "all",
smoother = TRUE,
qqbands = TRUE)
# `geom_smooth()` using formula 'y ~ x'
# `geom_smooth()` using formula 'y ~ x'
cooks1 <- cooks.distance(constraints.m1, group = "subj")
# Warning in cooks.distance.lmerMod(constraints.m1, group = "subj"): group is not
# a valid argument for this function. As of version 0.4.0, group has been replaced
# by level.
dotplot_diag(x = cooks1, cutoff = "internal", name = "cooks.distance", index=constraints$subj) + ylab("Cook's distance") + xlab("ID")
sjPlot::plot_model(constraints.m1,
type = "est",
colors = "bw",
sort.est=TRUE,
axis.title=c("Effect (Unexpected - Expected in seconds)"),
show.values=TRUE,
show.p=TRUE)
constraints.m1.table <- sjPlot:: tab_model(constraints.m1,
show.std =TRUE,
show.stat=TRUE,
show.df=TRUE)
# boundary (singular) fit: see ?isSingular
constraints.m1.beta <- summary(effectsize::standardize(constraints.m1))
# boundary (singular) fit: see ?isSingular
constraints.interventions <- cbind(
gen.beta(effectsize::standardize(constraints.m1)),
gen.m(constraints.m1),
gen.ci(constraints.m1)[3:10,]
)
# boundary (singular) fit: see ?isSingular
# Computing profile confidence intervals ...
constraints.cooks<- constraints[which(cooks1 <= 4/264),]
# where are the influential observations?
nrow(constraints)
# [1] 264
nrow(constraints.cooks)
# [1] 249
constraints.summary <- constraints %>%
group_by(paper,condition) %>%
summarise(n = n())
# `summarise()` regrouping output by 'paper' (override with `.groups` argument)
constraints.cooksn <- constraints.cooks %>%
group_by(paper,condition) %>%
summarise(n_cooks = n())
# `summarise()` regrouping output by 'paper' (override with `.groups` argument)
constraints.where.influential <- full_join(constraints.summary,constraints.cooksn)
# Joining, by = c("paper", "condition")
constraints.m1.cooks <- lmer(data = constraints.cooks,
formula = look_pref ~ training_yesno + action_causal + action_consequence + actor_hand + agent_efficient_fam + ageday + (1|condition) + (1|experiment) + (1|paper),
REML=FALSE)
# boundary (singular) fit: see ?isSingular
summary(constraints.m1.cooks)
# Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
# method [lmerModLmerTest]
# Formula: look_pref ~ training_yesno + action_causal + action_consequence +
# actor_hand + agent_efficient_fam + ageday + (1 | condition) +
# (1 | experiment) + (1 | paper)
# Data: constraints.cooks
#
# AIC BIC logLik deviance df.resid
# 1484.8 1527.0 -730.4 1460.8 237
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -2.9748 -0.6479 -0.0429 0.6569 3.1180
#
# Random effects:
# Groups Name Variance Std.Dev.
# condition (Intercept) 0.0 0.00
# experiment (Intercept) 0.0 0.00
# paper (Intercept) 0.0 0.00
# Residual 20.7 4.55
# Number of obs: 249, groups: condition, 12; experiment, 5; paper, 2
#
# Fixed effects:
# Estimate Std. Error df t value Pr(>|t|)
# (Intercept) -4.0078 3.6877 249.0000 -1.09 0.278
# training_yesno1 -1.9535 0.4773 249.0000 -4.09 0.00005773 ***
# action_causal1 -2.5634 0.4768 249.0000 -5.38 0.00000017 ***
# action_consequence1 -1.2255 0.6202 249.0000 -1.98 0.049 *
# actor_hand1 0.7581 0.8187 249.0000 0.93 0.355
# actor_hand2 0.9867 0.8311 249.0000 1.19 0.236
# agent_efficient_fam1 -1.9383 0.4261 249.0000 -4.55 0.00000844 ***
# ageday 0.0281 0.0333 249.0000 0.84 0.400
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Correlation of Fixed Effects:
# (Intr) trnn_1 actn_cs1 actn_cn1 actr_1 actr_2 agn__1
# tranng_ysn1 -0.061
# action_csl1 0.161 0.076
# actn_cnsqn1 0.007 0.058 0.330
# actor_hand1 0.120 -0.226 0.001 -0.377
# actor_hand2 -0.047 -0.223 0.002 0.744 -0.644
# agnt_ffcn_1 0.124 0.313 0.248 0.190 0.000 0.001
# ageday -0.983 0.018 -0.070 -0.035 -0.009 -0.025 -0.050
# convergence code: 0
# boundary (singular) fit: see ?isSingular
constraints.m1.cooks.beta <- lmer(data = constraints.cooks,
formula = scale(look_pref) ~ training_yesno + action_causal + action_consequence + actor_hand + agent_efficient_fam + scale(ageday) + (1|condition) + (1|experiment) + (1|paper),
REML=FALSE)
# boundary (singular) fit: see ?isSingular
summary(constraints.m1.cooks.beta)
# Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
# method [lmerModLmerTest]
# Formula:
# scale(look_pref) ~ training_yesno + action_causal + action_consequence +
# actor_hand + agent_efficient_fam + scale(ageday) + (1 | condition) +
# (1 | experiment) + (1 | paper)
# Data: constraints.cooks
#
# AIC BIC logLik deviance df.resid
# 683.5 725.7 -329.8 659.5 237
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -2.9748 -0.6479 -0.0429 0.6569 3.1180
#
# Random effects:
# Groups Name Variance Std.Dev.
# condition (Intercept) 0.000 0.00
# experiment (Intercept) 0.000 0.00
# paper (Intercept) 0.000 0.00
# Residual 0.828 0.91
# Number of obs: 249, groups: condition, 12; experiment, 5; paper, 2
#
# Fixed effects:
# Estimate Std. Error df t value Pr(>|t|)
# (Intercept) -0.3453 0.1338 249.0000 -2.58 0.010 *
# training_yesno1 -0.3908 0.0955 249.0000 -4.09 0.00005773 ***
# action_causal1 -0.5129 0.0954 249.0000 -5.38 0.00000017 ***
# action_consequence1 -0.2452 0.1241 249.0000 -1.98 0.049 *
# actor_hand1 0.1517 0.1638 249.0000 0.93 0.355
# actor_hand2 0.1974 0.1663 249.0000 1.19 0.236
# agent_efficient_fam1 -0.3878 0.0853 249.0000 -4.55 0.00000844 ***
# scale(ageday) 0.0490 0.0580 249.0000 0.84 0.400
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Correlation of Fixed Effects:
# (Intr) trnn_1 actn_cs1 actn_cn1 actr_1 actr_2 agn__1
# tranng_ysn1 -0.238
# action_csl1 0.510 0.076
# actn_cnsqn1 -0.151 0.058 0.330
# actor_hand1 0.613 -0.226 0.001 -0.377
# actor_hand2 -0.393 -0.223 0.002 0.744 -0.644
# agnt_ffcn_1 0.415 0.313 0.248 0.190 0.000 0.001
# scale(agdy) -0.058 0.018 -0.070 -0.035 -0.009 -0.025 -0.050
# convergence code: 0
# boundary (singular) fit: see ?isSingular
plot_model(constraints.m1.cooks.beta)
sjPlot::plot_model(constraints.m1.cooks,
type = "pred",
# colors = "bw",
sort.est=TRUE,
axis.title=c("Effect (Unexpected - Expected in seconds)"),
show.values=TRUE,
show.p=TRUE)
# $training_yesno
#
# $action_causal
#
# $action_consequence
#
# $actor_hand
#
# $agent_efficient_fam
#
# $ageday
plot(allEffects(constraints.m1.cooks))
sjPlot::plot_model(constraints.m1.cooks,
type = "est",
colors = "bw",
sort.est=TRUE,
axis.labels= rev(c("mittened", "gloved", "age", "state change", "actor initially rational", "sticky mittens", "action effect on contact")),
axis.title=c("Effect (Unexpected - Expected in seconds)"),
show.values=TRUE,
show.p=TRUE)
constraints.interventions.cooks <- cbind(
gen.beta(effectsize::standardize(constraints.m1.cooks)),
gen.m(constraints.m1.cooks),
gen.ci(constraints.m1.cooks)[3:10,]
)
# boundary (singular) fit: see ?isSingular
# Computing profile confidence intervals ...
When we compared the effects of different variants we found that motor training increased infants’ looking preference for the unexpected event ([4.175,4.978], ß=-0.391, B=-1.953, SE=0.477, p<.001, two-tailed). We also found that other interventions on the actions infants see, including seeing an action that causes an effect on contact ([-11.263,3.247], ß=-0.513, B=-2.563, SE=0.477, p<.001, two-tailed), and seeing a reaching action that results in an object changing state, rather than picking up that object ([-2.892,-1.014], ß=-0.245, B=-1.226, SE=0.62, p=0.049, two-tailed). These analyses included 249/264 total infants - the other 15 were classified as influential using Cook’s Distance. Including all subjects in the analysis produced similar results, except that the state change effect crossed the threshold into non-significance ([-3.382,-0.989], ß=-0.173, B=-1.069, SE=0.776, p=0.17, two-tailed). We did not find an effect of age, [-0.649,2.623], ß=0.049, B=0.028, SE=0.033, p=0.4, two-tailed.
goals <- ind.data %>% filter(task =="goals")
goals.baseline.data <- goals %>%
filter(condition %in% c("1_control",
"1_twoobject") &
paper %in% c("gerson2014a",
"luo2011"))
goals.b1 <- lmer(data = goals.baseline.data,
formula = look_pref ~ 1 + ageday + (1|condition))
goals.b1.table <- sjPlot:: tab_model(goals.b1,
show.std =TRUE,
show.stat=TRUE,
show.df=TRUE)
goals.b1.beta <- summary(effectsize::standardize(goals.b1))
goals.baseline <- cbind(
gen.beta(effectsize::standardize(goals.b1)),
gen.m(goals.b1),
gen.ci(goals.b1)[3:4,]
)
# Computing profile confidence intervals ...
cooks1 <- cooks.distance(goals.b1, group = "subj")
# Warning in cooks.distance.lmerMod(goals.b1, group = "subj"): group is not a
# valid argument for this function. As of version 0.4.0, group has been replaced
# by level.
dotplot_diag(x = cooks1, cutoff = "internal", name = "cooks.distance", index=goals.baseline.data$subj) + ylab("Cook's distance") + xlab("subjID")
excluded.subs <- c("2", "7", "6", "11")
goals.b1.cooks <- lmer(data = goals.baseline.data %>%
filter(!subj %in% excluded.subs),
formula = look_pref ~ 1 + ageday + (1|condition))
summary(goals.b1.cooks)
# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula: look_pref ~ 1 + ageday + (1 | condition)
# Data: goals.baseline.data %>% filter(!subj %in% excluded.subs)
#
# REML criterion at convergence: 239.2
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -1.6626 -0.5931 0.0995 0.4841 2.8202
#
# Random effects:
# Groups Name Variance Std.Dev.
# condition (Intercept) 202 14.2
# Residual 106 10.3
# Number of obs: 32, groups: condition, 2
#
# Fixed effects:
# Estimate Std. Error df t value Pr(>|t|)
# (Intercept) 44.720 27.063 21.547 1.65 0.11
# ageday -0.355 0.236 29.073 -1.51 0.14
#
# Correlation of Fixed Effects:
# (Intr)
# ageday -0.925
goals.baseline.cooks <- cbind(
gen.beta(effectsize::standardize(goals.b1.cooks)),
gen.m(goals.b1.cooks),
gen.ci(goals.b1.cooks)[3:4,]
)
# Computing profile confidence intervals ...
# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep
# Warning in profile.merMod(object, which = parm, signames = oldNames, ...): non-
# monotonic profile for (Intercept)
# Warning in confint.thpr(pp, level = level, zeta = zeta): bad spline fit for
# (Intercept): falling back to linear interpolation
For standard versions of the goals task, we found differential looking between the expected and unexpected events, Mean_diff = 3.918 seconds, [-0.931,-0.081], ß=-0.318, B=-0.503, SE=0.214, p=0.025, two-tailed. However, this was driven by 4 influential observations - without these 4 observations, there was no longer a reliable looking preference, [-0.815,0.126], ß=-0.21, B=-0.355, SE=0.236, p=0.143, two-tailed.
goals.m1 <- lmer(data = goals,
formula = look_pref ~ training_yesno + action_consequence + location_object_goal_ambiguous + agent + bothobjects_present_visible_fam + ageday + (1|condition) + (1|experiment) + (1|paper),
control = lmerControl(optimizer="Nelder_Mead"))
# boundary (singular) fit: see ?isSingular
summary(goals.m1)
# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula:
# look_pref ~ training_yesno + action_consequence + location_object_goal_ambiguous +
# agent + bothobjects_present_visible_fam + ageday + (1 | condition) +
# (1 | experiment) + (1 | paper)
# Data: goals
# Control: lmerControl(optimizer = "Nelder_Mead")
#
# REML criterion at convergence: 2842
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -4.049 -0.458 -0.028 0.416 2.921
#
# Random effects:
# Groups Name Variance Std.Dev.
# condition (Intercept) 19.6 4.43
# paper (Intercept) 0.0 0.00
# experiment (Intercept) 0.0 0.00
# Residual 153.4 12.39
# Number of obs: 362, groups: condition, 13; paper, 6; experiment, 2
#
# Fixed effects:
# Estimate Std. Error df t value Pr(>|t|)
# (Intercept) 10.3237 6.2887 96.2681 1.64 0.104
# training_yesno1 -1.9519 2.4419 4.8329 -0.80 0.462
# action_consequence1 2.8811 2.4462 7.9680 1.18 0.273
# location_object_goal_ambiguous1 5.1289 2.4944 8.5843 2.06 0.071
# agent1 -5.4390 2.9960 6.3026 -1.82 0.117
# agent2 6.4666 2.7646 19.3741 2.34 0.030
# bothobjects_present_visible_fam1 -5.3013 2.0586 15.0765 -2.58 0.021
# ageday -0.0733 0.0531 352.2521 -1.38 0.168
#
# (Intercept)
# training_yesno1
# action_consequence1
# location_object_goal_ambiguous1 .
# agent1
# agent2 *
# bothobjects_present_visible_fam1 *
# ageday
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Correlation of Fixed Effects:
# (Intr) trnn_1 actn_1 lc___1 agent1 agent2 bt___1
# tranng_ysn1 -0.295
# actn_cnsqn1 -0.246 0.002
# lctn_bjc__1 0.067 0.002 0.482
# agent1 0.011 0.547 -0.169 0.111
# agent2 -0.013 -0.297 0.024 -0.170 -0.796
# bthbjct___1 0.119 0.002 -0.245 -0.170 0.387 -0.563
# ageday -0.875 0.043 0.053 0.038 -0.042 0.017 0.057
# convergence code: 0
# boundary (singular) fit: see ?isSingular
goals %>%
select(paper, experiment, condition, training_yesno, object_diff_size_huge,action_consequence,location_object_goal_ambiguous,agent,bothobjects_present_visible_fam) %>%
distinct()
# paper experiment condition training_yesno
# 1 luo2011 1 1_twoobject no
# 2 luo2011 1 1_oneobject no
# 3 luo2011 2 2_differentpositions no
# 4 gerson2014a 1 1_active yes
# 5 gerson2014a 1 1_observation no
# 6 gerson2014a 1 1_control no
# 7 gerson2014b 1 1_active yes
# 8 gerson2014b 1 1_observation no
# 9 gerson2014b 2 2_generalization yes
# 10 choi_unpub 1 1_equidistant no
# 11 choi_unpub 1 1_neartarget no
# 12 choi_unpub 1 1_fartarget no
# 13 choi2018 1 1_twoobject no
# 14 choi2018 1 1_oneobject no
# 15 choi2018 1 1_hidden no
# 16 woo_unpub 1 1_statechange_objectswap no
# 17 woo_unpub 2 2_disambiguatingobjectgoal no
# object_diff_size_huge action_consequence location_object_goal_ambiguous
# 1 no none yes
# 2 no none yes
# 3 no none no
# 4 no none yes
# 5 no none yes
# 6 no none yes
# 7 no none yes
# 8 no none yes
# 9 no none yes
# 10 yes none yes
# 11 yes none yes
# 12 yes none yes
# 13 yes none yes
# 14 yes none yes
# 15 yes none yes
# 16 no state_change yes
# 17 no state_change no
# agent bothobjects_present_visible_fam
# 1 object yes
# 2 object no
# 3 object no
# 4 hand yes
# 5 hand yes
# 6 hand yes
# 7 hand yes
# 8 hand yes
# 9 hand yes
# 10 person yes
# 11 person yes
# 12 person no
# 13 person yes
# 14 person yes
# 15 person yes
# 16 person yes
# 17 person yes
# originally pre-registered model was rank deficient, needed to drop 2 variables
# dropped action_causal because it is redundant with action_consequence (i.e. there are no non-causal actions that are state changes)
# dropped object_diff_size_huge also because it is almost completely redundant with agent (i.e. all studies with a full human agent also are "yes" for "object_diff_size_huge")
# aa <- allFit(goals.m1)
summary(goals.m1)
# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula:
# look_pref ~ training_yesno + action_consequence + location_object_goal_ambiguous +
# agent + bothobjects_present_visible_fam + ageday + (1 | condition) +
# (1 | experiment) + (1 | paper)
# Data: goals
# Control: lmerControl(optimizer = "Nelder_Mead")
#
# REML criterion at convergence: 2842
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -4.049 -0.458 -0.028 0.416 2.921
#
# Random effects:
# Groups Name Variance Std.Dev.
# condition (Intercept) 19.6 4.43
# paper (Intercept) 0.0 0.00
# experiment (Intercept) 0.0 0.00
# Residual 153.4 12.39
# Number of obs: 362, groups: condition, 13; paper, 6; experiment, 2
#
# Fixed effects:
# Estimate Std. Error df t value Pr(>|t|)
# (Intercept) 10.3237 6.2887 96.2681 1.64 0.104
# training_yesno1 -1.9519 2.4419 4.8329 -0.80 0.462
# action_consequence1 2.8811 2.4462 7.9680 1.18 0.273
# location_object_goal_ambiguous1 5.1289 2.4944 8.5843 2.06 0.071
# agent1 -5.4390 2.9960 6.3026 -1.82 0.117
# agent2 6.4666 2.7646 19.3741 2.34 0.030
# bothobjects_present_visible_fam1 -5.3013 2.0586 15.0765 -2.58 0.021
# ageday -0.0733 0.0531 352.2521 -1.38 0.168
#
# (Intercept)
# training_yesno1
# action_consequence1
# location_object_goal_ambiguous1 .
# agent1
# agent2 *
# bothobjects_present_visible_fam1 *
# ageday
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Correlation of Fixed Effects:
# (Intr) trnn_1 actn_1 lc___1 agent1 agent2 bt___1
# tranng_ysn1 -0.295
# actn_cnsqn1 -0.246 0.002
# lctn_bjc__1 0.067 0.002 0.482
# agent1 0.011 0.547 -0.169 0.111
# agent2 -0.013 -0.297 0.024 -0.170 -0.796
# bthbjct___1 0.119 0.002 -0.245 -0.170 0.387 -0.563
# ageday -0.875 0.043 0.053 0.038 -0.042 0.017 0.057
# convergence code: 0
# boundary (singular) fit: see ?isSingular
sjPlot::plot_model(goals.m1,
type = "est",
colors = "bw",
sort.est=TRUE,
axis.title=c("Effect (Unexpected - Expected in seconds)"),
show.values=TRUE,
show.p=TRUE)
sjPlot:: tab_model(goals.m1, show.stat=TRUE)
| look pref | ||||
|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p |
| (Intercept) | 10.32 | -2.00 – 22.65 | 1.64 | 0.101 |
| training_yesno1 | -1.95 | -6.74 – 2.83 | -0.80 | 0.424 |
| action_consequence1 | 2.88 | -1.91 – 7.68 | 1.18 | 0.239 |
| location_object_goal_ambiguous1 | 5.13 | 0.24 – 10.02 | 2.06 | 0.040 |
| agent1 | -5.44 | -11.31 – 0.43 | -1.82 | 0.069 |
| agent2 | 6.47 | 1.05 – 11.89 | 2.34 | 0.019 |
| bothobjects_present_visible_fam1 | -5.30 | -9.34 – -1.27 | -2.58 | 0.010 |
| ageday | -0.07 | -0.18 – 0.03 | -1.38 | 0.168 |
| Random Effects | ||||
| σ2 | 153.43 | |||
| τ00 condition | 19.60 | |||
| τ00 paper | 0.00 | |||
| τ00 experiment | 0.00 | |||
| N condition | 13 | |||
| N experiment | 2 | |||
| N paper | 6 | |||
| Observations | 362 | |||
| Marginal R2 / Conditional R2 | 0.104 / NA | |||
plot(allEffects(goals.m1))
goals.interventions<- cbind(
gen.beta(effectsize::standardize(goals.m1)),
gen.m(goals.m1),
gen.ci(goals.m1)[4:11,]
)
# boundary (singular) fit: see ?isSingular
# Computing profile confidence intervals ...
# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep
# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep
# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep
# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep
# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep
# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep
# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep
# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep
# Warning in profile.merMod(object, which = parm, signames = oldNames, ...): non-
# monotonic profile for (Intercept)
# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep
# Warning in profile.merMod(object, which = parm, signames = oldNames, ...): non-
# monotonic profile for agent2
# Warning in confint.thpr(pp, level = level, zeta = zeta): bad spline fit for
# (Intercept): falling back to linear interpolation
# Warning in confint.thpr(pp, level = level, zeta = zeta): bad spline fit for
# agent2: falling back to linear interpolation
cooks1 <- cooks.distance(goals.m1, group = "subj")
# Warning in cooks.distance.lmerMod(goals.m1, group = "subj"): group is not a
# valid argument for this function. As of version 0.4.0, group has been replaced
# by level.
dotplot_diag(x = cooks1, cutoff = "internal", name = "cooks.distance", index=goals$subj) + ylab("Cook's distance") + xlab("subjID")
goals.cooks <- goals[which(cooks1 <= 4/362),]
# where are the influential observations?
nrow(goals)
# [1] 362
nrow(goals.cooks)
# [1] 345
goals.summary <- goals %>%
group_by(paper,condition) %>%
summarise(n = n())
# `summarise()` regrouping output by 'paper' (override with `.groups` argument)
goals.cooksn <- goals.cooks %>%
group_by(paper,condition) %>%
summarise(n_cooks = n())
# `summarise()` regrouping output by 'paper' (override with `.groups` argument)
goals.where.influential <- full_join(goals.summary, goals.cooksn)
# Joining, by = c("paper", "condition")
goals.m1.cooks <- lmer(data = goals.cooks,
formula = look_pref ~ training_yesno + action_consequence + location_object_goal_ambiguous + agent + bothobjects_present_visible_fam + ageday + (1|condition) + (1|experiment) + (1|paper),
control = lmerControl(optimizer="Nelder_Mead"))
# boundary (singular) fit: see ?isSingular
fixef(goals.m1.cooks)
# (Intercept) training_yesno1
# 6.06716 -1.83931
# action_consequence1 location_object_goal_ambiguous1
# 2.33221 5.15730
# agent1 agent2
# -4.46125 5.60015
# bothobjects_present_visible_fam1 ageday
# -3.07025 -0.01669
fixef(goals.m1.cooks,add.dropped=TRUE)
# (Intercept) training_yesno1
# 6.06716 -1.83931
# action_consequence1 location_object_goal_ambiguous1
# 2.33221 5.15730
# agent1 agent2
# -4.46125 5.60015
# bothobjects_present_visible_fam1 ageday
# -3.07025 -0.01669
summary(goals.m1.cooks)
# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula:
# look_pref ~ training_yesno + action_consequence + location_object_goal_ambiguous +
# agent + bothobjects_present_visible_fam + ageday + (1 | condition) +
# (1 | experiment) + (1 | paper)
# Data: goals.cooks
# Control: lmerControl(optimizer = "Nelder_Mead")
#
# REML criterion at convergence: 2595
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -3.484 -0.459 -0.042 0.468 3.371
#
# Random effects:
# Groups Name Variance Std.Dev.
# condition (Intercept) 13.3 3.64
# paper (Intercept) 0.0 0.00
# experiment (Intercept) 0.0 0.00
# Residual 110.0 10.49
# Number of obs: 345, groups: condition, 13; paper, 6; experiment, 2
#
# Fixed effects:
# Estimate Std. Error df t value Pr(>|t|)
# (Intercept) 6.0672 5.4877 109.6470 1.11 0.271
# training_yesno1 -1.8393 2.0190 4.7758 -0.91 0.406
# action_consequence1 2.3322 2.0883 8.2716 1.12 0.295
# location_object_goal_ambiguous1 5.1573 2.2202 9.8410 2.32 0.043
# agent1 -4.4613 2.5736 6.7760 -1.73 0.128
# agent2 5.6002 2.6245 23.8723 2.13 0.043
# bothobjects_present_visible_fam1 -3.0702 1.8673 12.4200 -1.64 0.125
# ageday -0.0167 0.0465 336.0909 -0.36 0.720
#
# (Intercept)
# training_yesno1
# action_consequence1
# location_object_goal_ambiguous1 *
# agent1
# agent2 *
# bothobjects_present_visible_fam1
# ageday
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Correlation of Fixed Effects:
# (Intr) trnn_1 actn_1 lc___1 agent1 agent2 bt___1
# tranng_ysn1 -0.284
# actn_cnsqn1 -0.226 0.002
# lctn_bjc__1 0.067 0.002 0.519
# agent1 -0.005 0.527 -0.154 0.123
# agent2 0.027 -0.260 -0.007 -0.201 -0.796
# bthbjct___1 0.114 0.003 -0.247 -0.151 0.411 -0.545
# ageday -0.880 0.046 0.048 0.034 -0.023 -0.014 0.073
# convergence code: 0
# boundary (singular) fit: see ?isSingular
plot(allEffects(goals.m1.cooks))
goals.interventions.cooks<- cbind(
gen.beta(effectsize::standardize(goals.m1.cooks)),
gen.m(goals.m1.cooks),
gen.ci(goals.m1.cooks)[4:11,]
)
# boundary (singular) fit: see ?isSingular
# Computing profile confidence intervals ...
We note that the analysis we reported deviates from our pre-registration plan. In that plan, we included 8 fixed effects, including a fixed effect for age in days, one picking out a control condition, and 6 others. However, 2 of these fixed effects were either completely or partially redundant with other predictors in the model, which caused a rank deficiency issue. Thus, we dropped these two fixed effects from the model. We report this new model below.
When we compared the effects of different variants, we found that making the goal of the agent’s actions unambiguous ([-1.191,5.201], ß=0.465, B=5.157, SE=2.22, p=0.043, two-tailed), or presenting a self-propelled object as an agent ([-8.285,-0.729], ß=0.505, B=5.6, SE=2.625, p=0.043, two-tailed), relative to a hand or a person, increased infants’ looking preference for the unexpected event. We did not find any other effects, including those of motor experience, or action consequences, when taking into account all of these predictors in the same model. These analyses included 345/362 total infants - the other 17 were classified as influential using Cook’s Distance. Including all subjects in the analysis produced similar results, except that the goal ambiguity effect crossed the threshold into marginal-significance ([-3.382,-0.989], ß=-0.173, B=-1.069, SE=0.776, p=0.17, two-tailed). We did not find an effect of age, [-5.763,0.015], ß=-0.019, B=-0.017, SE=0.047, p=0.72, two-tailed.